2022 대한지리학회 연례학술대회 포스터 세션

국지 온도를 고려한 서울특별시 공공 무더위 쉼터 접근성 분석 - 고령층을 중심으로

※ 본 문서는 PC 환경에 최적화되어있습니다.

서울대학교 지리학과 SGIS 학회에서 진행한 프로젝트 결과물로, 대한지리학회 포스터 세션에서 발표한 연구이다.

1 서론

제6차 IPCC 평가 보고서에 따르면, 최근 2,000년 동안 전례 없는 속도의 지표 온도 상승이 나타나고 있다. 2021년 한국의 평균기온은 1973년 이후 두 번째로 높은 평균기온 13.3℃를 기록했다. 이러한 온도 상승은 인류의 삶과 생태계에 큰 위협으로 작용하며, 특히 폭염은 태풍, 호우, 폭설 등의 자연 재난과 달리 그 피해가 주로 건강 및 생명과 직결된다는 점에서 위협적이다.
개인의 특성에 따라 폭염의 영향을 받는 정도는 상이하며, 폭염에 취약한 대표적인 계층으로 고령층을 꼽을 수 있다. 당뇨병, 심혈 질환, 호흡기 질환, 신장 질환 등의 기저질환을 앓는 고령층은 특히 폭염에 더 취약한데, 36.5℃ 이상의 환경에서 기온이 1℃씩 상승할 때마다 65세 이상 인구의 사망률은 최대 28.4%까지 증가한다. 한국은 세계에서 고령화가 급속도로 진행되고 있는 나라 중 하나로, 고령화 현상은 폭염 재해 취약계층의 증가를 의미하며, 이들 집단을 위한 기후 대비의 필요성이 제기된다.
이러한 폭염 재해는 토지 개발 패턴, 열에 대한 물리적 노출 정도, 사회적 취약성 등이 지역마다 이질적으로 분포한다는 점에서 공간성에 대한 고려가 중요하게 작용한다. 예를 들어, 이상 기후와 관련된 연구에 있어 지표 기온 측정에 일반적으로 활용되는 원격탐사는 광범위한 지역의 자료를 손쉽게 제공하여 현상의 분석에 유용한 시각을 가져다주는 이점이 있다. 이러한 지리공간적 기법의 응용은 위험의 시공간적 변화 양상에 대한 깊은 이해를 가능케 한다.
그러나 폭염 재해의 위험성에 관한 연구는 국내에서의 연구 진행이 미흡한 실정이다. 따라서 본 연구에서는 고령 인구의 공간적 분포 양상, 서울특별시의 국지적 기온 분포 양상, 무더위 쉼터의 운영 현황을 각각 살펴보고, 이들을 종합하여 고령층을 중심으로 한 서울특별시의 폭염 재해 취약지역을 도출하고자 한다. 서울특별시는 이번 5월부터 2022 여름철 종합대책을 추진하며, 올여름 무더위쉼터 3,400여 곳을 운영한다고 발표한 바 있다. 이러한 사회적 분위기에 발맞춰, 본 연구를 통해 서울특별시의 폭염 재해 정책에 대한 의사결정에 시사점을 제안해보고자 한다.

2 방법론

본 연구에서는 취약성 평가를 위하여 열 환경과 폭염 민감성이라는 두 변수를 설정하고, 이변량 국지모란지수를 통해 두 변수의 공간적 군집을 탐색, 고령층을 중심으로 한 서울특별시 폭염 취약 지역을 산출한다.

2.1 개념 설정

고령층을 대상으로 한 폭염 취약 지역의 산출에 앞서, 먼저 ’폭염 취약성’의 개념에 대한 정의가 선행되어야 한다. 취약성은 인구, 사회 시스템 또는 특정 장소가 위험 및 재해에 얼마나 취약한지를 의미하는데, 이와 관련한 다수의 선행 연구에서 폭염 취약성은 폭염에의 노출 정도, 폭염에 대한 반응성, 폭염에 대한 적응 가능성으로 측정된 바 있다. 폭염에의 노출 정도와 폭염에 대한 반응성은 폭염 취약성을 유발하는 인자이며, 폭염에 대한 적응 가능성은 폭염 취약성을 감쇄하는 인자라고 할 수 있다. 본 연구에서 폭염에의 노출 정도는 지표면 온도를 기반으로 한 열 환경으로, 폭염에 대한 반응성은 고령인구의 수로, 폭염에 대한 적응 가능성은 공공 무더위쉼터의 이용 가능성으로 평가한다.

2.2 지표 선정

2.2.1 열 환경

지표면 온도(LST)는 보건 및 열 취약성 연구에 널리 사용되어 왔으며, 본 연구에 사용된 LST 데이터는 100m 공간해상도를 지니는 TRIS 감지기(Thermal Infrared Sensor)로 획득한 Landsat-8 Collection 2 Level 2 영상에서 도출되었다. 미 지질조사국(USGS)은 이미 원본 영상의 밝기값을 Landsat Collection 2 Level-1의 열적외선 대역, 대기상단(TOA) 반사율, TOA 밝기 온도, ASTER(Advanced Spaceborne Thermal Emission and Reflection Radiometer) GED(Global Emissivity Database) 데이터, ASTER 정규식생지수(NDVI) 데이터, Goddard Earth Observing System(GEOS)에서 추출한 지오퍼텐셜 높이, 비습 및 기온 등을 사용하여 변환한 지표온도 자료를 제공하고 있으며 본 연구는 이를 사용하였다. 또한 구한 접근성 지수와 스케일을 동일하게 하기 위하여 행정동 경계 안에 위치한 셀을 평균화하는, 즉 동 평균화 작업을 한 자료를 사용하였다. 시점은 온열질환이 주로 발생하는 여름철(6~8월)을 잡았으며 지표면 온도 추정에 심각한 영향을 미칠 수 있는 구름이 10% 이하인 영상으로 분석을 실시하였다.

2.2.2 폭염 민감성

다음으로, 폭염 민감성에 대한 변수의 경우 그 지표는 2SFCA(Two-step Floating Catchment Area)기반의 접근성 지수를 이용한다. Radke and Mu에 의해 처음 제안되어 Luo and Wang에 의해 수정된 2SFCA는 중력모델 기반의 공간적 접근성 측정 방법으로, 공급시설에 대한 수요자 대비 공급능력과 수요 임계거리 내 서비스 공급시설의 인구대비 공급자 비율을 합산하여 접근성을 산출하며 다음의 수식으로 표현할 수 있다.

\[ R_j=\frac{S_j}{\sum_{k\in\{d_{kj}\le d_0\}}D_kf(d_{kj})} \]

\[ A_i=\sum_{j\in\{{d_{ij}\le d_0\}}}R_jf(d_{ij}) \\ A_i= \sum_j\frac{S_jf(d_{ij})}{\sum_{k\in\{d_{kj}\le d_0\}}D_kf(d_{kj})} \]

여기서 \(A_i\)는 \(i\)지점에서의 접근성이고,\(R_j\)는 \(j\)지점에 위친한 공급시설의 공급-수요 비, \(S_j\)는 \(j\)지점의 공급시설 가중치(혹은 공급량),\(D_k\)는 \(k\)지점의 수요, \(d_0\)는 이동비용 한계거리 \(d_{ij}\) 혹은 \(d_{kj}\)은 각각 \(i\),\(j\) 혹은 \(k\),\(j\) 간의 거리, \(f()\)는 거리조락함수를 의미한다.

본 연구에서는 서울특별시 공공 무더위쉼터의 위치정보와 면적, 행정동별 65세 고령인구 수를 활용하여 65세 이상 고령인구의 수를 수요로, 공공 무더위쉼터의 면적을 공급으로, 노인 보행 특성을 고려한 300m를 수요 임계거리로 설정한 뒤 분석을 진행하였다.

2.3 공간적 패턴 및 군집 탐색

폭염 취약성과 관련한 두 변수 사이의 관계와 공간적 패턴 및 군집 양상을 살펴보기 위해 이변량 국지모란지수를 활용하였다. Anselin et al.에 의해 개발된 이변량 국지모란지수(Bivariate Local Indicators of Spatial Association)는 서로 다른 종류의 변수들 사이에서 특정 관측치 주변에 유사한 값들이 형성하는 공간적 집적의 유의도를 측정하는 기법으로, 두 변수 간의 국지적 공간 상관성을 측정하고 공간 상관관계 유형을 식별하는데 활용된다. 이는 다음과 같이 정의된다.

\[ I_{kl}=Z_k^i\sum_{j=1}^nW_{ij}Z_l^j \]

여기에서 \(I_{kl}\)은 이변량 국지모란 값을 의미하며, \(Z_{ik}\)는 기준 지역 \(I\)에서 변수 \(k\)의 값을, \(Z_{jl}\)은 기준 지역 \(j\)에서 변수 \(l\)의 값을 나타낸다. \(W_{ij}\)는 \(i\)지역과 인접한 \(j\)지역에 관한 공간가중행렬을 의미한다. 이를 통해 기준 지역 \(i\)에서 변수 \(k\)의 값과 기준 지역과 인접한 \(j\)지역의 변수 \(l\)의 값을 비교할 수 있다. 이 값은 양(음)의 값은 양(음)의 상관관계가 있음을 의미하거나 값이 유사한 관측치의 집중(분산)이 있음을 의미한다. 그 결과로 4가지 패턴을 추출할 수 있는데, 기준 지역의 값이 높고 주변에 이와 비슷한 값이 집중해있음을 나타내는 High-High, 기준 지역의 값이 높고 주변에 이와 다른 값이 집중해있음을 나타내는 High-Low, 기준 지역의 값이 낮고 주변에 이와 비슷한 값이 집중해있음 나타내는 Low-Low, 기준 지역의 값이 낮고 주변에 이와 다른 값이 집중해 있음을 나타내는 Low-High 패턴이 존재한다. 본 연구에서는 고령인구의 무더위쉼터 접근성이 낮은 곳이 폭염에의 취약 정도가 높다는 측면을 강조하고자 2SFCA 접근성 지수에 음수를 취하였다. 이를 통해 접근성이 낮은, 즉 취약도가 높은 곳에 인접한 지역에서 높은 지표면 온도가 나타나면 HIgh-High 패턴이 나타날 수 있도록 하였다.

일련의 연구 과정을 도식화하면 아래와 같다.

3 데이터

3.1 위성영상 자료

USGS에서 제공하는 Landsat-8 위성 원격탐사 자료를 활용한 LST 표면을 사용하였다. 이 때 관측일시가 2019년 6월 13일인 원격탐사자료를 활용하였는데 이는 여름철 달인 6월~8월 중 하나를 충족하며, 지표구름피복 정도가 10% 미만이면서 OLI 및 TIRS 센서의 영상자료 질이 9(-1: 품질 미평가, 0: 최악, 9: 최고)로서 분석에 용이하다고 판단하였기 때문이다. 이에 본 조는 해당 자료를 분석 시 자료로 활용하였다.

3.2 서울시 무더위쉼터 현황

서울시 무더위쉼터 현황은 서울 열린데이터광장에서 제공하는 ‘서울시 무더위쉼터’ 데이터를 활용하였다. 2021년 하절기에 운영된 경로당, 사회복지관, 주민센터의 공공 무더위쉼터를 대상으로 하였으며, 총 3,049개의 무더위쉼터에 대한 이름, 도로명주소, 위치정보, 면적을 구득하였다.

3.3 65세 이상 고령인구 현황

고령층, 즉 65세 이상 고령인구 현황은 서울 열린데이터광장에서 제공하는 서울시 고령자 현황 (동별) 통계 데이터를 활용하였다. 2021년 3/4분기 기준의 동별 65세 이상 고령자 수를 구득하였다.

3.4 데이터 정제 과정

3.4.1 서울특별시 shp 로드 및 정제

library(sf)
seoulemd <- read_sf("./data/shp/Z_SOP_BND_ADM_DONG_PG.shp",options = "ENCODING=euc-kr")[c(1:425),] %>% 
  st_make_valid() %>%
  rmapshaper::ms_simplify(keep = 0.05,
                          keep_shapes = FALSE) %>%
  st_make_valid() %>%
  st_transform(32652) %>%
  st_cast("MULTIPOLYGON")

3.4.2 고령인구 데이터 로드 및 병합

library(tidyverse)
pop <- read_tsv("./data/demo/report.txt", skip = 2) %>%
  filter(동 != "합계" & 동 != "소계") %>% 
  dplyr::select("기간", "자치구", "동", "65세이상고령자")


names(pop)[4] <- c("고령인구수")
pop[is.na(pop)] <- 0
pop$동 <- gsub("\\.","·", pop$동)
pop_gangil <- rbind(pop[424, 4], pop[426, 4]) # 강일동 = 강일동 + 상일제2동, 상일동 = 상일제1동
pop_gangil <- colSums(pop_gangil)
pop <- rbind(pop[1:423, ], c("2021.3/4", "강동구", "강일동", pop_gangil), c("2021.3/4", "강동구", "상일동", 4326))




code <- read_csv("./data/demo/adm_code.csv", skip = 1) %>% filter(시도코드 == 11) 
code[, c(2,3,4,5)] <- code[, c(3,5,2,4)]
code <- code %>% 
  unite(행정동코드, 시도코드, 시도명칭, 시군구코드) %>%
  dplyr::select(-시군구명칭)
code$행정동코드 <- gsub("_","", code$행정동코드)
colnames(code) = c("ADM_DR_CD", "자치구", "동")


pop <- pop %>% 
  full_join(code, by = c("자치구", "동"))
pop[, 4] <- pop %>% 
  dplyr::select("고령인구수") %>%  
  mutate_all("as.numeric")


seoulemd <- seoulemd %>% 
  full_join(pop, by = "ADM_DR_CD") %>% 
  dplyr::select("ADM_DR_CD", "기간", "자치구", "동", "고령인구수")

3.4.3 무더위쉼터 데이터 로드 및 정제

target <- c(1,2,5)

hotshelter <- read_csv("./data/shelter/hotshelter.csv") %>%
  filter(시설유형 %in% target) %>%
  st_as_sf(coords = c("경도", "위도"),
           crs = 4326) %>%
  st_transform(st_crs(seoulemd))

3.4.4 시각화

# 무더위쉼터 매핑
library(tmap)

tm_shape(seoulemd) + tm_polygons() + 
 tm_scale_bar(position = c("left", "bottom")) +
 tm_shape(hotshelter) + tm_dots(size = 0.01) +
 tm_layout(main.title = "서울특별시 무더위쉼터 분포 현황", main.title.size = 1.2, frame = FALSE)

# 동별 고령인구 단계구분도
tm_shape(seoulemd) + tm_polygons(col = "고령인구수") +
 tm_scale_bar(position = c("left", "bottom")) +
 tm_layout(main.title = "서울특별시 동별 고령인구 수", main.title.size = 1.2, frame = FALSE)

3.4.5 고령인구의 무더위쉼터 공간 접근성

3.4.5.1 Points in Polygon

seoulhs <- hotshelter %>% 
  st_join(seoulemd) %>%
  group_by(동) %>% 
  summarize(hotshelter = n())


seoulhs <- st_drop_geometry(seoulhs)

seoulemd <- seoulemd %>%
  left_join(seoulhs, by = "동") %>%
  mutate(hotshelter = replace_na(hotshelter, 0))

seoulemd <- seoulemd %>%
  mutate(hsperpop = (hotshelter/고령인구수) * 1000)


tm_shape(seoulemd) + tm_polygons(col = "hsperpop", style = "jenks", palette = "PuBu", border.alpha = 0)

3.4.5.2 Point Proximity Buffers

centroid_seoulemd <- st_centroid(seoulemd)

ct_buff <- st_buffer(centroid_seoulemd, dist = 300)

hs_buff <- hotshelter %>% st_join(ct_buff) %>% group_by(동) %>% summarize(hs300m = n())
hs_buff <- st_drop_geometry(hs_buff)
seoulemd <- seoulemd %>% 
 left_join(hs_buff, by = "동") %>% 
 mutate(hs300m = replace_na(hs300m, 0)) %>% 
 mutate(hsbuff300m = (hs300m/고령인구수) * 1000)

ggplot(seoulemd) + geom_histogram(mapping = aes(x = hs300m), na.rm = TRUE)

tm_shape(seoulemd) + tm_polygons(col = "hsbuff300m", style = "jenks", palette = "PuBu", border.alpha = 0)

3.4.5.3 Distance to Nearest Site

library(matrixStats)

hotshelter_dist <- st_distance(centroid_seoulemd, hotshelter)
seoulemd <- seoulemd %>% mutate(hotsheltermin = rowMins(hotshelter_dist))

tm_shape(seoulemd) +
  tm_polygons(col = "hotsheltermin", style = "jenks", palette = "PuBu", border.alpha = 0)

3.4.5.4 Floating Catchment Area (FCA)

library(SpatialAcc)

buff_tract <- centroid_seoulemd %>% 
        dplyr::select(고령인구수) %>% 
        st_join(ct_buff) %>%
        group_by(동) %>% 
        summarize(buffpop = sum(고령인구수.x)) %>%
        ungroup()
buff_tract <- st_drop_geometry(buff_tract)

seoulemd <- seoulemd %>%
                      left_join(buff_tract, by = "동") %>%
                      mutate(buffpop = replace_na(buffpop, 0)) %>%
                      mutate(FCA = (hs300m/buffpop) * 10000)

tm_shape(seoulemd) +
  tm_polygons(col = "FCA", style = "jenks", palette = "PuBu", border.alpha = 0)

3.4.5.5 Two-step Floating Catchment Area (2SFCA)

centroid_coords <- st_coordinates(centroid_seoulemd)
hotshelter_coords <- st_coordinates(hotshelter)
dist_matrix <- SpatialAcc::distance(centroid_coords, hotshelter_coords, type = "manhattan")

TSFCA <- ac(p = seoulemd$고령인구수, 
          n = hotshelter$면적,
          D = dist_matrix, d0 = 300, family = "2SFCA")
seoulemd <- seoulemd %>% mutate(TSFCA = TSFCA)
seoulemd <- seoulemd %>% mutate(`-TSFCA` = (-1*TSFCA))

tm_shape(seoulemd) +
 tm_polygons(col = "TSFCA", palette = "RdBu", style = "jenks", border.alpha = 0, legend.hist = TRUE) +
 tm_layout(legend.outside = TRUE)

3.4.6 서울시 LST 표면 생성

library(stars)

B10 <- read_stars("./data/raster/landsat/LC08_L2SP_116034_20190613_20200828_02_T1_ST_B10.TIF")


# 참고
# https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/media/files/LSDS-1619_Landsat-8-9-C2-L2-ScienceProductGuide-v4.pdf

LST <- (B10*0.00341802) + 149 - 273.15


names(LST) = "Land Surface Temperature of Seoul - 20190613 1110hrs in deg C"


plot(LST[seoulemd], breaks = "equal", col = hcl.colors(11, rev = T, "Spectral"))

3.4.7 데이터 표면 동평균화

위 과정에서 구한 여러 데이터 표면을 행정동 경계를 사용하여 동 평균을 한 뒤, 이를 행정동 자료의 속성정보로 넣어주는 과정이다.

3.4.8 토지피복도 불러오기

서울시 중앙을 가로지르는 한강은 동 평균시에 근처의 행정동에 막대한 영향을 끼쳐 한강을 비롯한 수역을 제거하였다.

다음과 같은 코드로 R 상에서 세분류 토지피복도의 모든 서울 도엽을 불러와 작업하였다. 많은 셰이프 파일을 R로 불러오기 위해 병렬컴퓨팅을 사용하였다.

#```{r eval=FALSE, include=TRUE}

## Land Cover Classification
fi <-  list.files(path = "./data/shp/토지피복/", pattern = "*.shp", full.names = T)



lcccrs <- read_sf(fi[1]) %>% st_crs()



library(doParallel)  

{
no_cores <- detectCores() - 1  
cl <- makeCluster(no_cores, type="PSOCK")  
registerDoParallel(cl)  
LCC <- foreach(i=1:length(fi),
               #.combine = rbind,
               .packages = c("sf")) %dopar% {
                 read_sf(fi[i])
               }

LCC <- mapedit:::combine_list_of_sf(LCC) %>% st_transform(st_crs(seoulemd))


stopCluster(cl)
}


#st_write(LCC,"./test.shp",append=F)

3.4.9 수역 제외 동평균 과정

수역만 제외하고자 L1 코드가 700인 데이터만 추출하였다.

#수계만 추출
water <- filter(LCC, L1_CODE == "700")

masked <- raster::mask(as(LST[seoulemd],"Raster") , water, inverse = TRUE)


raster::plot(masked)

다음은 수역을 제외한 나머지 기온 데이터표면을 동평균 자료로 바꿔주는 작업이다.

#수계를 제외한 온도표면 산출
library(velox) # 고속 래스터 연산 패키지
library(raster)

lst.land <- velox(masked)


# 평균
m.tmp <- lst.land$extract(as(seoulemd,"Spatial"), fun = function(x) mean(x, na.rm = TRUE))

#행정동 자료에 넣기
pcsseoul <-  cbind(seoulemd,m.tmp)


tmap_mode("plot")
tm_shape(pcsseoul) + tm_polygons(title = "19.06.14. 1110시 서울시 행정동별 평균지표온", 
              id="sggnm", 
              col = "m.tmp", 
              #popup.vars = c("adm_cd","sggnm","adm_nm","m.tmp"),
              style = "jenks",
              n = 7,
              palette = "-RdBu",
              legend.hist = TRUE
  ) +
  tm_layout(legend.outside = T)

#library(ggplot2)
#g1 = ggplot() + geom_sf(data = pcsseoul, aes(fill = m.tmp)) 
#g1

3.4.10 이변량 모란

R 을 이용하여 일련의 분석을 진행한 결과, 강서·양천구 일대, 금천구 일대, 동대문·성동·광진구 일대, 강동·송파구 일대에서 통계적으로 유의한 High-High 패턴 클러스터 4개를 식별할 수 있었다. 이들 클러스터는 고령층의 공공 무더위 쉼터 접근 취약성이 높은 동시에 인근 지표면 온도 역시 높은 곳으로, 해당 지역 에 우선적인 공공 무더위 쉼터의 확충이 필요함을 시사한다.

이변량 모란 계산 방법은 spdep 패키지를 활용한 사용자 정의 함수를 만들어 계산하는 방법과 Rgeoda 패키지를 사용한 방법이 있다.

3.4.10.1 User-defined Function

library(spdep)
library(spatialreg)
library(Matrix)
#======================================================
# Programming some functions

# Bivariate Moran's I
moran_I <- function(x, y = NULL, W){
        if(is.null(y)) y = x

        xp <- (x - mean(x, na.rm=T))/sd(x, na.rm=T)
        yp <- (y - mean(y, na.rm=T))/sd(y, na.rm=T)
        W[which(is.na(W))] <- 0
        n <- nrow(W)

        global <- (xp%*%W%*%yp)/(n - 1)
        local  <- (xp*W%*%yp)

        list(global = global, local  = as.numeric(local))
}


# Permutations for the Bivariate Moran's I
simula_moran <- function(x, y = NULL, W, nsims = 1000){

        if(is.null(y)) y = x

        n   = nrow(W)
        IDs = 1:n

        xp <- (x - mean(x, na.rm=T))/sd(x, na.rm=T)
        W[which(is.na(W))] <- 0

        global_sims = NULL
        local_sims  = matrix(NA, nrow = n, ncol=nsims)

        ID_sample = sample(IDs, size = n*nsims, replace = T)

        y_s = y[ID_sample]
        y_s = matrix(y_s, nrow = n, ncol = nsims)
        y_s <- (y_s - apply(y_s, 1, mean))/apply(y_s, 1, sd)

        global_sims  <- as.numeric( (xp%*%W%*%y_s)/(n - 1) )
        local_sims  <- (xp*W%*%y_s)

        list(global_sims = global_sims,
             local_sims  = local_sims)
}


#======================================================
# Adjacency Matrix (Queen)

nb <- poly2nb(pcsseoul, queen=T, snap=2)
lw <- nb2listw(nb, style = "B", zero.policy = T)
W  <- as(lw, "symmetricMatrix")
W  <- as.matrix(W/rowSums(W))
W[which(is.na(W))] <- 0

#======================================================
# Calculating the index and its simulated distribution
# for global and local values

m   <- moran_I(pcsseoul$X.TSFCA, pcsseoul$m.tmp, W)
m[[1]] # global value

            [,1]
[1,] -0.08836747

m_i <- m[[2]]  # local values

local_sims <- simula_moran(pcsseoul$X.TSFCA, pcsseoul$m.tmp, W)$local_sims

# Identifying the significant values 
alpha <- .05  # for a 95% confidence interval
probs <- c(alpha/2, 1-alpha/2)
intervals <- t( apply(local_sims, 1, function(x) quantile(x, probs=probs)))
sig        <- ( m_i < intervals[,1] )  | ( m_i > intervals[,2] )

#======================================================
# Preparing for plotting

#oregon.tract     <- st_as_sf(oregon.tract)
pcsseoul$sig <- sig


x <- pcsseoul$X.TSFCA
y <- pcsseoul$m.tmp

# Identifying the LISA patterns
xp <- (x-mean(x))/sd(x)
yp <- (y-mean(y))/sd(y)

patterns <- as.character( interaction(xp > 0, W%*%yp > 0) ) 
patterns <- patterns %>% 
        str_replace_all("TRUE","High") %>% 
        str_replace_all("FALSE","Low")
patterns[pcsseoul$sig==0] <- "Not significant"
pcsseoul$patterns <- patterns




# Plotting
ggplot() + geom_sf(data=pcsseoul, aes(fill=patterns), color="NA") +
        scale_fill_manual(values = c("red", "pink", "light blue", "dark blue", "grey95")) + 
        theme_minimal()

library(tmap)
tmap_mode("view")
  
tm_shape(pcsseoul) + tm_polygons(title = c("Bivar LISA","hsperpop","TSFCA","Mean Temperatre"), 
              id="동", 
              col = c("patterns","hsperpop","TSFCA","m.tmp"),
              popup.vars = c("자치구","동","고령인구수","hsperpop","hs300m","TSFCA","m.tmp","patterns"),
              style = "jenks",
              alpha = 0.6
              )+
  tm_facets(sync = TRUE, ncol = 2)

3.4.10.1.1 Rgeoda

#혹은 Rgeoda
library(rgeoda)

wij <- queen_weights(precision_threshold = 2, pcsseoul)

( qsa <- local_bimoran(wij, pcsseoul[c('X.TSFCA', 'm.tmp')]) )

Reference class object of class "LISA"
Field "gda_lisa":
An object of class "p_LISA"
Slot "pointer":
<pointer: 0x0000026f476e4be0>

Field "p_vals":
 [1] 0.238 0.175 0.001 0.001 0.035 0.211 0.258 0.182 0.001 0.047 0.001 0.003
[13] 0.002 0.002 0.001 0.011 0.352 0.454 0.282 0.302 0.048 0.062 0.013 0.157
[25] 0.001 0.001 0.308 0.001 0.147 0.205 0.040 0.004 0.274 0.231 0.094 0.400
[37] 0.117 0.272 0.169 0.363 0.328 0.364 0.470 0.396 0.110 0.006 0.155 0.379
[49] 0.173 0.039 0.207 0.151 0.464 0.353 0.447 0.290 0.408 0.476 0.267 0.187
[61] 0.007 0.009 0.001 0.341 0.243 0.017 0.110 0.013 0.245 0.147 0.195 0.046
[73] 0.218 0.446 0.485
 [ reached getOption("max.print") -- omitted 350 entries ]
Field "c_vals":
 [1] 0 0 4 4 2 0 0 0 3 1 3 1 3 3 3 2 0 0 0 0 1 0 3 0 3 3 0 1 0 0 1 3 0 0 0 0 0 0
[39] 0 0 0 0 0 0 0 3 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 1 0 1 0 0 0 1 0 0 0
 [ reached getOption("max.print") -- omitted 350 entries ]
Field "lisa_vals":
 [1] -0.165584762 -0.215592579 -0.274538500 -1.068847536  2.878158356
 [6] -0.033014326 -0.197450421 -0.195419049 -0.485906130  0.401741651
[11] -0.640319932  0.115604118 -6.778421814 -0.556526797 -1.565842005
[16]  5.393412726 -0.054049316 -0.004296734 -0.458209217  0.131993063
[21]  0.339921846  0.303438401 -0.940771218 -0.490937833 -0.309274783
[26] -0.827225870 -0.008206477  0.694310939 -0.034610060 -0.001892288
[31]  0.006683266 -0.469768512  0.185798031  0.140990861 -0.183322714
[36] -0.228935596 -0.072989938 -2.580975541 -0.201118383 -0.085117112
[41] -0.554581446 -0.140299427  0.021021477 -0.023299325  0.267280108
[46] -0.619766515  0.056706520  0.068672254 -0.250413710  0.316821351
[51]  0.198745471 -0.988637185 -0.004296273 -0.023099662 -0.002786939
[56]  0.004361840  0.045096883  0.015727250  0.013093975  0.017268080
[61]  0.451959132  0.492454208  0.459327477 -0.038526296 -0.119201772
[66]  0.385972958  0.265265719  0.156439509 -0.015344593 -0.588438437
[71]  0.191530294  0.332508201 -0.004154887 -0.009060773  0.003789772
 [ reached getOption("max.print") -- omitted 350 entries ]
Field "nn_vals":
 [1] 6 6 6 7 6 5 3 9 7 5 7 5 4 5 5 7 6 6 7 6 7 7 5 5 5 6 5 9 6 5 5 7 4 8 6 6 5 4
[39] 7 5 5 4 9 4 6 6 7 9 6 5 6 6 4 6 6 6 7 6 4 6 8 7 6 6 7 8 4 6 4 5 6 7 6 5 6
 [ reached getOption("max.print") -- omitted 350 entries ]
Field "labels":
[1] "Not significant" "High-High"       "Low-Low"         "Low-High"       
[5] "High-Low"        "Undefined"       "Isolated"       
Field "colors":
[1] "#eeeeee" "#FF0000" "#0000FF" "#a7adf9" "#f4ada8" "#464646" "#999999"

  lisa_colors <- lisa_colors(qsa)
  lisa_labels <- lisa_labels(qsa)
  lisa_clusters <- lisa_clusters(qsa)

{  plot(st_geometry(pcsseoul), col=sapply(lisa_clusters, 
     function(x){return(lisa_colors[[x+1]])}), 
     border = "#333333", lwd=0.2)
title(main = "Bivaraite Local Moran")
legend('bottomleft', legend = lisa_labels, fill = lisa_colors, 
       border = "#eeeeee")
  }

library(tmap)
tmap_mode("view")
  
tm_shape(pcsseoul) + tm_polygons(title = c("Bivar LISA","hsperpop","TSFCA","Mean Temperatre"), 
              id="동", 
              col = c("patterns","hsperpop","TSFCA","m.tmp"),
              popup.vars = c("자치구","동","고령인구수","hsperpop","hs300m","TSFCA","m.tmp","patterns"),
              style = "jenks",
              alpha = 0.6
              )+
  tm_facets(sync = TRUE, ncol = 2)

outseoul <- pcsseoul[,c(1:4)]

outseoul$Pattern <- lisa_clusters
outseoul$pval <- lisa_pvalues(qsa)


for (i in 1:length(outseoul$ADM_DR_CD)) {
  j <- as.numeric(outseoul$Pattern[i])+1
  outseoul$Pattern[i] <- lisa_labels[j]
}


ggplot() + geom_sf(data=outseoul, aes(fill=Pattern), color="NA") +
        scale_fill_manual(values = c("red", "pink", "light blue", "dark blue", "grey95")) + 
        theme_minimal()+ ggtitle("서울특별시 공공 무더위 쉼터 접근성 분석 결과")+theme(plot.title=element_text(face="bold",size=20))

library("classInt")
library("RColorBrewer")

brks<-classIntervals(outseoul$pval, n=5, style="fixed", fixedBreaks=c(0.01, 0.03, 0.05, 0.1, 1)) #define categories
brks <- round(brks$brks,digits=2) #round
catVar<-findInterval(outseoul$pval, brks, all.inside=TRUE) #assign categories for the data (previously melted using 'melt' function)

outseoul$`p-value` <- catVar #join data & spatial info

# Create labels from break values
intLabels <- matrix(1:(length(brks)-1))
for(i in 1:length(intLabels )){intLabels [i] <- paste(as.character(brks[i]),"-",as.character(brks[i+1]))}
#intLabels[1,1]<-c("< 0.8")
intLabels[4,1]<-c(">0.1")

#actual plotting step
ggplot() + geom_sf(data=outseoul, aes(fill=`p-value`), color="NA") +
 scale_fill_gradientn(colours=rev(brewer.pal(4, "RdBu")), guide="legend", label=intLabels) + theme_minimal()+ggtitle("이변량 국지모란 유의수준 분포")+theme(plot.title=element_text(face="bold",size=20))

4 결론 및 토의

최근 심화되는 이상기후와 기온 상승 문제는 한국 사회의 고령화와 맞물려 폭염재해가 고령층에게 큰 위험으로 작용하는 배경이 되었다. 이에, 본 연구에서는 서울특별시를 중심으로 고령층의 폭염재해에 대한 취약성을 평가하고, 취약 지역을 도출하였다. 본 연구의 한계 및 개선 방안은 다음과 같다. 먼저, 수요에 해당하는 고령인구의 구체적인 분포를 고려하지 못하였다. 대시메트릭 매핑 등의 기법을 이용하여 고령인구의 인구 중심점을 더욱 현실에 맞게 수정한다면 더욱 정확한 연구를 진행할 수 있을 것이다. 둘째, 자료 접근의 한계로 인하여 폭염에 대한 반응성과 적응 가능성을 평가하는 변수로 한정적인 지표만을 활용하였다. 행정동 단위의 소득 수준, 냉방기기 보급률을 평가할 수 있는 자료가 확보된다면 추후 이를 활용하여 연구의 정확성이 높아질 수 있으리라 기대한다.